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Abstract 

We propose and analyze an effective free energy describing the 
physics of dischnation defects in particle arrays constrained to move 
on an arbitrary two-dimensional surface. At finite temperature the 
physics of interacting disclinations is mapped to a Laplacian Sine- 
Gordon Hamiltonian suitable for numerical simulations. We then spe- 
cialize to the case of a spherical crystal at zero temperature. The 
ground state is analyzed as a function of the ratio of the defect core 
energy to the Young's modulus. We argue that the core energy con- 
tribution becomes less and less important in the limit R ^ a, where 
R is the radius of the sphere and a is the particle spacing. For large 
core energies there are twelve disclinations forming an icosahedron. 
For intermediate core energies unusual finite-length grain boundaries 
are preferred. The complicated regime of small core energies, appro- 
priate to the limit R/a oo, is also addressed. Finally we discuss the 
application of our results to the classic Thomson problem of finding 
the ground state of electrons distributed on a two-sphere. 
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1 Introduction 



The theory of two-dimensional melting of essentially planar materials (mono- 
layers) is a rich and well-developed subject [jl], ^f. An interesting aspect of 
melting in this low dimension is that both the crystalline to hexatic and hex- 
atic to fluid transitions can be driven by the sequential liberation of point-like 
topological defects - dislocations in the former case and disclinations in the 
latter. It is clearly important, therefore, to have a thorough understanding 
of the statistical mechanics of interacting topological defects. On the plane 
all topological defects are bound at zero temperature, but on manifolds with 
more complicated topology excess free disclinations must exist even at zero 
temperature. 

The statistical mechanics of particles confined to frozen surfaces of con- 
stant positive and negative curvature was discussed, e.g. in references 
and [Q. It was argued that regions of positive and negative curvature would 
promote the formation of unpaired disclinations, and that these might be 
screened by clouds of dislocations. At low temperature, it was suggested 
that the anisotropic interaction between these screening dislocations would 
lead them to condense into grain boundaries. The physics of particles on a 
quenched random topography was discussed in Ref. [Q. 

The simplest example of a surface with positive Gaussian curvature is the 
sphere. Dodgson studied the ground state of the Abrikosov flux lattice in a 
model thin film superconductor on a sphere (subject to a field radiating from 
a magnetic monopole at the center), and found evidence for twelve five- fold 
disclination defects at the vertices of an icosahedron in an otherwise six- 
coordinated crystalline environment 0. This defect configuration is similar 
to one proposed by Lubensky and collaborators for lipid bilayer vesicles in the 
hexatic phase 0, except that in hexatics the disclination energy is reduced 
by screening due to an equilibrium concentration of unbound dislocations. 
Later, Dodgson and Moore proposed adding dislocations to the ground state 
of a sufficiently large vortex crystal in a spherical geometry to screen out the 
strains associated with twelve extra disclinations in the Abrikosov phase . 
Vortices in a thin film superconductor behave like particles interacting with 
a repulsive logarithmic pair potential. Another context in which crystalline 
ground states on a sphere arise is the so-called Thomson problem, where the 
vortices are replaced by particles interacting with a repulsive 1/r potential 



10, 11]. Our own interest in this class of problems was stimulated by 



the beautiful work of Alar Toomre W2\, which we discuss later (and which 
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hopefully will be described by Toomre himself one day!). Toomre's ideas also 
play a key role in a recent paper on the Thomson problem by Perez-Garrido 
and Moore |T^. For a discussion of disclination and dislocation defects for 
disk-like configurations of electrons in the plane see [Q. 

The study of melting and the nature of the ground state on curved man- 
ifolds may be a good testing ground for understanding the new features that 
arise from the topological defects required for particle arrays on surfaces with 
nontrivial topologies. Our approach is to work directly with the defects them- 
selves, and treat the particles within continuum elastic theory. This approach 
is more general than, say, a direct simulation of particles interacting with a 
logarithmic or 1/r potential, because all details of the pair potential are em- 
bodied in the elastic constants mediating the interaction between defects. 
By eliminating explicit reference to the particles themselves, we also greatly 
reduce the number of degrees of freedom needed to study the ground state. 
As we shall see, the effective Hamiltonian used here, in which defects such 
as grain boundaries and dislocations are built up out of elementary disclina- 
tions, leads to a variety of interesting and novel structures not encountered 
in the plane. 

The statistical mechanics of monolayers on curved surfaces such as the 
sphere may also be viewed as the infinite bending rigidity limit of membranes 
with a spherical topology. Our investigation may therefore be considered a 
prelude to the careful incorporation of defects in the study of the phase 
transitions of, e.g., membranes composed of lipid bilayers | |15|| . 

It is useful to review expectations for low temperature configurations of 
crystals in flat space |jl6l . Although the ground state is believed to be defect 



free, one can certainly consider the response to adding a single excess discli- 
nation. The stresses induced by such a disclination are very high, and the 
energy can be lowered by polarizing the surrounding medium into dislocation 
pairs, as indicated schematically in Fig. |l[ 

When interactions between dislocations are taken into account one might 
expect them to organize into grain boundaries (i.e., lines of dislocations with 
Burgers vectors oriented perpendicular to the lines) to minimize the en- 
ergy even further. Experiments on smectic liquid crystal films with tilted 
molecules (the tilt is used to force in an extra disclination) reveal a pat- 
tern of five jagged grain boundaries radiating outward, consistent with this 
picture. Computer simulations with periodic boundary conditions have been 
used to study the relaxation of a disclination quartet (two fives and two sev- 
ens), from an initial configuration where these defects sit on the corners of a 
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Figure 1: Schematic of an isolated threefold disclination, in an approximately 
four-coordinated medium. The elastic stress in the vicinity of the isolated 
disclination is relieved by the formation of a screening cloud of dislocations. 



very large square in an otherwise six-coordinated medium |T8[. After the re- 
laxation, grain boundaries joining the fives to the sevens appear. Relaxation 
of the disclination elastic stresses in this way occurs at a price - the core 
energies associated with the extra dislocations lead to an additional term in 
the energy which diverges linearly with system size R, as compared to the 
i?^ divergence associated with an unscreened disclination |1T9| . 

A situation reminiscent of these flat space experiments occurs on surfaces 
of non-zero Gaussian curvature, e.g., the sphere. Although the Gaussian 
curvature of the sphere approximately compensates the strains associated 
with isolated disclinations, for small core energies (or, equivalently, large 
sphere radius compared to the particle spacing) it can still be favorable to 
introduce extra screening dislocations into the ground state. 

To see how screening of an isolated five-fold disclination by dislocations 
comes out on a sphere, it is helpful to first consider what happens in fiat 
space. A five-fold disclination can be created by removing a wedge of material 
subtending an angle s = 2tt/6 and then deforming the remaining material to 
close the gap. (The disclination in the square lattice of Fig.|I| was made by 
removing a 27r/4 = 90° wedge). The resulting stresses were calculated, e.g., 
in Ref. pO[- We use polar coordinates r and 0, measured from the center 
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of the disclination. If /i and A are the material elastic constants, the stress 
tensor is dominated by a^^, where (neglecting logarithmic corrections due to 
boundary effects), 



and Kq, the Young's modulus, is related to the Lame coefficients by 

° 2/i + A ^ ^ 

Note that a^p^ is proportional to the disclination charge s. This approximately 
constant stress arises from the stretching of material required to close the 
gap engendered by the missing wedge, and leads to the divergence in 
disclination energy with system size R Consider now the fate of a 

tightly bound dislocation pair (the Burgers' vectors are equal and opposite 
with |6| = a) placed in the stress field of this disclination. The stress a^if, 



creates a Peach- Kohler force which tries to tear the dislocation pair apart |]2T 
We assume for simplicity a purely radial separation Ar between dislocations 
with Burgers' vectors in the tangential direction. The energy of the pair then 
consists of 2Ed {Ed is the dislocation core energy), a logarithmic binding 
energy and a linear Peach-Kohler term (proportional to a^rt>) which tries to 
"ionize" the pair, similar to the effect of an electric field on a charge dipole, 

Epa^r{Ar) = 2Ed + ^ln(— ) - ^s(Ar) . (3) 
The energy can be lowered once Ar exceeds Ar* ^ b/s ^ a and the pair 



separates |^ . One of the liberated dislocations moves off to infinity while the 
other remains to help screen the disclination. As more and more dislocations 
are created in this way, the stress is reduced until the dislocation density 



in an annulus of width dr at radius r from the disclination is ||2T[ 



^dW~(7^)— • (4) 
2tt ra 

Note that if these dislocations collapse to form a single linear grain boundary 
radiating out from the disclination, the angular deficit s is related to the 
spacing I between dislocations in the grain by s ~ a//. If the dislocations 
form m grain boundaries, the spacing will be / ~ am/s. In this paper we 
shall study the cases m = 2 and m = 5. 
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The analysis sketched above is easily adapted to the curved surface of a 
sphere of radius R. The angular deficit s associated with a circuit around a 
dislocation in flat space can now be compensated by the nonzero Gaussian 
curvature 1/R^. Let us assume that a five-fold disclination is placed at 
the north pole of the sphere. We describe the physics by geodesic polar 
coordinates (r, 0) about this point with metric 



ds' = dr^ + R\m'^{^)d(l)'^ , (5) 
R 

and work in the limit R^ a. We expect that the stress is controlled by the 
effective disclination charge inside a circuit at fixed geodesic distance r from 
the disclination (see Eq.(|l3]) below), namely 



^effir) = s-y^ d(l) dr'y/gK 



Note that Seff{r) decreases with increasing r. In the limit of weak curvature 
i? 3> a, we expect that the stress formula (|1]) is replaced by 

a,,{r) = Ko'-^ , (7) 

with a corresponding weakening of the Peach-Kohler force. The reduction in 
the angular stress o",^^ with increasing geodesic distance from the disclination 
arises because the stretching required to remove a wedge in fiat space is 
reduced according to the metric (^. We now expect the dislocations in m 
grain boundaries radiating from a disclination to exhibit a variable spacing 
between dislocations, 

Note that the spacing diverges as r — r ~ , where 

5 = e, = cos"^- = 33.56°. (9) 

The angular jump As(r) = Seff/m across the grain boundaries thus becomes 
smaller with increasing r and these boundaries eventually terminate when 
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the dislocations become sufficiently dilute. The total number of dislocations 
contained in the m grain boundaries is approximately 

Nd ~ ^s(x) . (10) 

Our calculations support this picture, and we find that the extra dislocations 
seem to form grain boundaries. Remarkably, and in contrast to flat space, 
these grain boundaries do indeed stop or start inside the crystalline medium. 
Our results also hint at a branching pattern of grain boundary networks (each 
radiating from a disclination), reminiscent of those found in Ref. ||17|| . 

The paper is organized as follows. In section we develop a formalism 
whose basic degrees of freedom are the defects themselves, rather than the 
underlying interacting particles. The particles themselves are treated within 
continuum elastic theory. As mentioned above, such a formalism has the ad- 
vantage of reducing the number of degrees of freedom as well as being rather 
universal in the sense that it applies to a broad class of interaction poten- 
tials. Varying the pair potential simply corresponds to changing the elastic 
moduli and defect core energy of the model. The model has the advantage 
that defects can move directly to positions which minimize the energy, with- 
out the constraints associated with disclination motion or dislocation climb 
in a crystalline medium which would attend a particle simulation. Despite 
its simplicity, finite temperature statistical mechanics of this model is still 
not amenable to a direct analytic solution. A duality mapping to an equiva- 
lent Laplacian Sine-Gordon model, however, yields a model with short range 
interactions whose lattice version should be straightforward to simulate nu- 
merically. 

In contrast, the limit of zero temperature may be treated analytically 
and we turn to this in the next three sections of the paper. In particular we 
discuss the ground state of a spherical crystal as a function of defect core 
energy relative to the combination of elastic constants (Young's modulus) 
which determines defect interactions at large distances. 

In section ^ we shift our attention from defects alone to the underlying 
lattice structure. We first discuss lattices with icosahedral symmetry. Our 
formalism applied to this case predicts the range of core energies for which 
the lattice is unstable to the formation of defects. 

An interesting application for our formalism is to the Thomson problem 
0, |T0[, discussed in section |^. The predictions of our approach are in agree- 



ment with existing results where comparisons are available. A beautiful ex- 
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perimental realization of the Thomson problem is provided by multi-electron 
bubbles trapped in liquid helium at low temperatures . Order in electrons 



confined by a positively charged capacitor plate to a helium surface has been 
studied for many years. Except for capillary wave deformations, crystalliza- 
tion proceeds in an essentially fiat environment. At high electron densities, 
curvature is introduced via an instability to a regular array of "dimples" in 
the helium surface, each containing a million electrons or more. Upon in- 
creasing the density of positive charge below the surface further by adding a 
metallic tip to the anode, one can form completely submerged multi-electron 
spherical bubbles. Typical bubbles contain 10^—10^ electrons. The outward 
electrostatic repulsion of the electrons on the inner surface of the helium 
bubble balances against the surface tension of the helium interface to pro- 
duce bubbles with diameters in the range 10—100 microns. Results for the 
Thomson problem have implications for trapped multi-electron bubbles well 
below the fiat space freezing temperature. 



2 Finite Temperature 
2.1 Free Energy 

As our main interest lies in the study of defects on two-dimensional curved 
surfaces, we need a formalism that deals directly with the defect degrees 
of freedom themselves. A rigorous geometrical derivation of the effective 
free energy for the defects is given in [24|. An equivalent derivation may 



also be given by integrating out the phonon degrees of freedom from the 
elastic Hamiltonian 



25 



E 



(fx 



with the appropriate modifications for a general 
distribution of defects. The energy of a two-dimensional crystal embedded 
in an arbitrary frozen geometry described by a metric giji'x) is given by 

1 

A2 
1 

A xy 

where g{x.) is the determinant of the metric tensor, K{x.) is the associated 
Gaussian curvature and s(x) the disclination density 



s(x)) 
- ^(x)) 



{K{y) - s{y)) 

r 

(K(y)-.(y))(ll) 



.(x) 
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with disclinations located at the sites Xj of an underlying triangulated 
particle array. The "charges" qi may be positive or negative. Although we 
do not restrict the allowed values of the charge, we expect the unit charge 
defects to dominate for energetic reasons. A plus one charge corresponds to a 
five-fold coordinated particle (a five-disclination) and a minus one charge cor- 
responds to a seven- fold coordinated particle (a seven-disclination). Charges 
are attracted to regions of like-sign Gaussian curvature. 

The first term of Eq. (|l2) represents a long range elastic interaction and 
Kq is the Young's modulus of Eq. (H) 

The second term in Eq. (|lT|) contains a single inverse-Laplacian operator, 
which is singular at short distances due to distortions of the lattice at dis- 
tances less than the lattice spacing. This is the dominant term for hexatic 



membranes, where Ka is the hexatic stiffness, [pq] In the present context, this 
singular contribution leads to a renormalized core energy Ecore{KA) for each 
defect and it represents non-universal details of the interaction on the scale 



of the inter-particle spacing a. The energy of Eq. (O) is thus simplified to 



+ NE,ore . (13) 

Although it is not essential, we assume for convenience that the core energies 
of five and seven-fold disclinations are identical. The partition function of 
our model is then 

m = E -ifwfy"''''"'- / n d^va n ^x^v/^e-^^^^^) , (m) 

where E{Kq) is the first term in Eq. (|TB|), y is the disclination fugacity e'^^"""" 
{(3 is the inverse temperature), N± is the total number of fives and sevens 
respectively and x is the Euler characteristic of the surface. For a given 
microscopic interaction potential both y and Kq are fixed. We shall find it 
useful, however, to regard these as independent parameters and discuss, in 
particular, the limits of large and small Ecore compared to Kqo^, where a is 
the lattice constant. 

Despite its elegant form this model is difficult to solve analytically. It is, 
moreover, challenging for direct numerical simulation because of the long- 
range interaction embodied in - see the explicit form for E{Kq) given 
in Eqs. (^8]) and (^9]) below. An alternative formulation is suggested by the 
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Laplacian roughening model for flat space melting . Direct molecular 

dynamics simulations or energy evaluations of particles interacting with a 
specified potential [12| are also of considerable interest. Since this ap- 
proach takes the particles as the primary degrees of freedom, rather than the 
defects, it falls outside the scope of the present paper. 



2.2 The Sine-Gordon model 

We now restrict ourselves to the case of the sphere, which has Euler char- 
acteristic X = 2. We map the previous model to a dual Sine-Gordon model 
with only short range interactions by adapting the derivation presented in 
to the present case, with some additional improvements. 
Let us start with the identity 

g-| / ^du^dv{s{u)-K{u))^{siv)~K{v)) 



The topological constraints ensure that the zero mode does not contribute 
to the path integral and this is indicated by the primes in the determinant 
and the measure. 

Since the zero mode is the constant eigenvector of the Laplacian, or- 
thonormality implies that 

J du^(t> = (16) 

for any configuration included in the measure of the path integral in 
Eq. (0). 

With this identity Eq. (|^) becomes 

Z{(3) = J p0'^(0)e-i/'^"v^^*^*e-*/'^"^'^^(") = J I?(/.'e-^W , (17) 

where the last identification defines the Hamiltonian 7f(0) and J^(0) is given 
by 



N+ N- 
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J: -^^^^(y J du^ef^^^-YHy J du^e^^^^-Y- (18) 
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Upon writing the Kronecker delta as 

Sn,-n^,u= C dx e-^'-<^--^-''\ (19) 

one finds ^ 

= [ dx e-*i2xg/d«V9cos(f0{«)+x) /2o) 
Jo 

Inserting this resuh in Eq. (^) and performing the integral over the x variable 
leads to 

n{ct>) = ^il?J du^A(PA<p-2y J du^ cos <P + ^ J du^Kct> . (21) 

The last (imaginary) term is a nuisance for practical applications. For the 
case of the sphere, however, the Gaussian curvature is constant, and we have 

j du^K(p = K J duy/g(p = , (22) 

where we have used Eq. (0). The Sine-Gordon representation for the sphere 
takes a very simple form 

W) = ^(^)'/ du^A(j)A<p-2y J duy^ cos(P . (23) 

Discretizing this expression for large y will yield a simple model with integer 
variables (f){u). A numerical simulation of this model seems the appropriate 
way to study the finite temperature statistical mechanics of defect arrays on 
a sphere ||26|, ^ . 

We now turn to the limiting case of zero temperature. 



3 Zero Temperature Limit 
3.1 General Surfaces 

The zero temperature limit requires the determination of the ground state 
by a minimization of the energy as a function of both the position and total 
number of defects. 

For the minimization with respect to the location of defects we see that 



the energy Eq. (13) depends only on the difference between the geometric 
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curvature and the defect density. As a result the defects will arrange them- 
selves to approximately match the Gaussian curvature determined by the 
geometry of the confining surface. A complete screening of the Gaussian 
curvature would yield a crystal with zero elastic energy at zero temperature. 
An important example is that of a crystal with the symmetry of a perfect 
icosahedron. The twelve positive disclinations located at its twelve vertices 
compensate the Gaussian curvature. There are twelve five-fold coordinated 
particles at the vertices, and all the rest are six-fold coordinated. 

As for the minimization with respect to total defect number, it is clear 
that the second term of Eq. (0) is linear with the number of defects, and 
so will clearly favor the lowest possible number of them. The physics of the 
zero temperature limit is therefore controlled by the competition between 
the core energy cost of creating a defect and the compensating gain from the 
screening of Gaussian curvature when defects are allowed to proliferate. 

3.2 The spherical crystal 

From now on we concentrate on a spherical crystal. Since the sphere has 
Euler characteristic 2 (genus 0) the charges qi of a set of disclinations must 
satisfy 

^ N 

J rf2xy^7(x)s(x) = 47r ^ ^ = 12 . (24) 

i=l 

This implies that, even at zero temperature, a sphere contains at least twelve 
excess five-fold disclinations. 

To evaluate the free energy Eq. (|TT|) we compute first the inverse square- 
Laplacian operator on a sphere of radius R, 

^x{OA,eA)-^-R E 1.^ ^^T^Yp ' ^^^^ 

where Yjj^{6, (p) are the spherical harmonics and {6, 0) are the usual spherical 
angles. The / = term does not appear in the sum, again as a result of 
the precise topology of the sphere (Eq. (^4|) ). The absence of this zero mode 
leads to a finite sum. The expression Eq. ( pS]) may also be written 

CO 9/ _|_ 1 

x(r, r; e\ <p') ^ x{f3) = R'j: ^,^^l ^^. Pii^osp) , (26) 
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where 

cos /3 = cos cos + sin 9" sin 9^ cos(0" - cf)^) (27) 

gives the length j3 of the geodesic arc connecting (^", 0") and (6'*, ip^) on the 
sphere. It is shown in Appendix P that this last sum may be written 



1 — cos/3 



In z 



(28) 



In Appendix ^ we discuss the flat space limit of infinite sphere radius. In 
Fig. 1^ we plot x/R^ (Eq. (^1)) as a function of the geodesic distance /3. 
Although the formula Eq. 



^ is simple, it is not particularly suitable for 
rapid numerical evaluation. In Appendix y we give alternative expressions 
for X better suited to fast numerical evaluation. 

The final expression for the total energy of a spherical crystal with an 
arbitrary number of disclinations follows from Eq. (B3) and Eq. fUT 



TX 



Kn 



N N 



36 



R'J2T.mxiO\<l)';9\<l)^) + NE, 



=1 j=i 



(29) 




Figure 2: Plot of x/R^ as a function of the geodesic angle /?. Only the 
interval /? G [0, vr] is plotted. 

Our interpretation of the disclination density screening out the Gaussian 
curvature can be made more precise. Note that 

TV 



.(x) 



'^V9 i=i 
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1 oo i TV 

= ;^ + ^E E >l(^><^)E9^>l*(^-<^^) (30) 

oo Z N 

= ^(x) + ^E E >l(^'</')E?^>t(^-</'0' 

1=1 m=—l i=l 

where the topological constraint Eq. (^) has been used. This last identity 
makes it clear that the set of equations 

AT 

5:g,F^(^„0,,)=O, (31) 

i=l 

for all / > 1 and all m, is the condition that the disclination density ex- 
actly matches the Gaussian curvature. Because it is difficult to imagine how 
discrete disclination charges could exactly cancel a smooth background Gaus- 
sian curvature, we cannot expect that Eq . (|3l|) will be satisfied in general for 
all values of /. We will, however, give examples where this set of equations 
is partially satisfied. It is easy to see, in fact, that in the limit of vanishing 
core energies, a configuration of defects satisfying Eq. (^Tj) is an absolute 
minimum of the energy Eq. (pQ]) , since the latter can be rewritten as 



^2 TV- oo ; Iv^ n-Y^ 

7C Kq x-^ \l^i=l'di^rr: 



E = E E " +ne^.^ . (32) 



2 



9 ^1^=.^ 

Eq. (^) then implies, for Ecore = 0, that the energy attains its minimum 
value of 0. 

Finally, we note that an equivalent expression for Eq. (p9D is given by 

ttKo.^ 21 + 1 ^ ^ 



E = ^ 12(1 E E mjPiicos A,) + N E^r, . (33) 

6K) L If + J-J i=i j=i 



Eqs. (|32D and (^91) are useful because they express the total energy as a sum 
of individual /-mode contributions 



oo 



E = Y.Ei, (34) 



1=1 



with the order of magnitude of each /-mode coefficient being roughly 
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By considering increasingly exotic arrangements of defects, we might hope 
to satisfy Eq. (^) for more and more low /-modes. If we do not enhance the 
large-Z contributions and do not pay too large a price in defect core energy, 
then the total energy will be small. 



4 Large core energies: The icosahedral lat- 
tice 

In the limit of large core energies the creation of additional defects will be 
strongly penalized and the sphere will contain only the minimum allowed 
twelve positive disclinations. From symmetry considerations it is a good 
ansatz to assume that these twelve disclinations minimize the repulsive x 
potential acting between them by forming an icosahedron X ||^, 0]. It is not 
difficult to check that the icosahedron is in fact an extremum of the energy 
Eq. m 



dE 



dE 
, -z- 



, (36) 



where i = 1, ■ ■ ■ , 12. We have checked numerically that fluctuations around 
this extremum increase the energy. Allowing the fluctuations to relax results 
in fast convergence to the icosahedron. Our numerical minimization gives the 
icosahedron as a global minimum. Thus our model successfully predicts an 
icosahedron minimum in the case where just twelve disclinations are allowed. 
From Eq. (|32D the energy is a function of the quantity 

Vi(I) = f:Yi{d,,<P,) , m = -l,---,l (37) 

i=l 

where (^i, 0i)i=i,...,i2 are particular coordinates for an icosahedron X on the 
sphere. This solution is obviously invariant under the full icosahedral group 
plus inversions, yh = yxCi. Since y is contained in SO (3), we can construct 
a representation of 3^ out of the irreducible representations of 5*0(3). We 
have 

Y: Dl^,{gy)Vi, = l-V:^, (38) 

m'=—l 

where gy is any element belonging to y. That is, \^ is a singlet of the 
icosahedral group 3^. This in turn means that if the trivial representation 
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(the so-called A representation) of the icosahedral group is not contained 
as an induced representation from the full rotational group, then Eq.(pID is 
exactly satisfied for any twelve disclinations forming an icosahedron, 

Vi = Y.Y^ie,,<P,) = 0, m = -l,---,l. (39) 

i 

It remains to identify those values of I which contain the trivial representa- 
tion. This is easily answered from an analysis of the characters of the group. 
The number of trivial representations ha contained in the /— th representa- 
tion of SO (3) is given by 

if sin{(/ + i)f} sin{(/ + i)f} 

60 I sm(|) sm(^) 

sin|(/ + i)f| r 1 .1 

20^-^ + 15sm{(/ + -).} , (40) 



which is nonzero for / = 6, 10, 12, 16 and all even-/, / > 16 |31[]. Note that 
Eq. (^) is satisfied for all /-odd modes, as follows from applying the inversion 
operator /, the generator of the Ci subgroup of 3^h- 

The icosahedral solution screens out the Gaussian curvature very effec- 
tively. Eq. (^) is partially satisfied, particularly for low /. The icosahedral 
lattice allows for non-zero contributions for only three (/ = 6, 10, 12) of the 
first fifteen putative contributions in Eq. (|30|) . A numerical evaluation gives 
the energy of an icosahedron E-^ as 

= 0.604 i^R') + 12 E,,re , (41) 
36 

where E^ore is the core energy of a single disclination. Its precise value is 
non-universal and depends on short distance details of the microscopic pair 
potential. The coefficient of ttKqB? is universal, independent of short- 
distance properties. Let us study it in more detail. In Table |l| we show the 
relative contribution from each /-mode. It is apparent that the first allowed 
non-zero contribution / = 6 accounts for almost 80% of the total energy of 
the icosahedron. Note also the relatively rapid convergence of the expansion; 
truncating up to the / = 100 mode gives a result which differs by less than 
0.2% from the actual result, Eq. (|4l|). It is remarkable how much the energy 
would be reduced by canceling out the / = 6 mode without further enhancing 
higher /-modes. 
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/ Ef Ef I Ef Ef 



6 


0.4669 


0.4669 


30 


0.0017 


0.5925 


10 


0.0329 


0.4999 


40 


0.0012 


0.5975 


12 


0.0507 


0.5506 


50 


0.0004 


0.5997 


16 


0.0129 


0.5635 


80 


0.0001 


0.6025 


18 


0.0125 


0.5760 


100 


3 X 10-^ 


0.6031 


20 


0.0004 


0.5764 


oo 





0.6043 



Table 1: The first column is the particular mode considered. Ef is the 
contribution of the /-mode to the total energy and Ef the running sum after 
adding all modes less than or equal to /. For convenience we set Ecore = 
and 7tKqR^/3Q = 1 in this table. 



5 Small core energies: The proliferation of 
defects 

If the defect core energies are small then the elastic energy may be reduced 
by creating additional defects. The topological constraint Eq. ( ^41) requires 
that additional defects appear in pairs of opposite charge. The challenge 
now is to understand and study the different possible distributions of these 
charges and the reduction in energy that those bring about when compared 
with the pure icosahedral case. The general form of the energy is, similarly 
to Eq. (grp, 



^ = C— ^ - +NE,ore, (42) 

36 V a / 

where we introduce the C-coefficient as a convenient parametrization of the 
elastic part of the energy. In the limit R/a — >■ oo, where a is the par- 
ticle spacing, we expect grain boundaries containing N ~ R/a disloca- 
tions emerging from each disclination. Hence the elastic term will always 
dominate over the core energy term in this limit. The critical sphere ra- 
dius R = Rc above which long range elastic energies dominate is given by 

Rc ~ const. (36-E'coT-e/7r-ft'oC^^)0'- 

If the total number of defects is large, an unconstrained minimization of 
Eq. (p9|) becomes an involved numerical problem. Instead of pursuing this 
further, we develop different approximations that allow us to tackle the case 
of a large number of defects while still capturing the most important features 
of the problem. 
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With results for a particular defect ansatz expressed as in Eq. (^^D, we 
need to determine the minimum distance of closest approach of neighboring 
plus-minus defect pairs. Upon identifying this distance with a, we can obtain 
the total number of particles M embodied in the defect configuration via the 
identification 

M«-|jfl^ (43) 
5.1 The icosahedral approximation 

Let us add new sets of twelve defects, each set lying on the vertices of an 
icosahedron. That is, we consider Eq. (^91) , not as a function of individual 
defects, but as a function of icosahedra of defects. From the mathematical 
arguments in the previous section, we can guarantee that the /-modes which 
vanish in the expansion Eq. ( P^ for the pure icosahedral case, will continue 
to do so within this approach. Since most of the low Z-modes, which dominate 
the energy, vanish for any icosahedron, we expect that the Euler angles of 
the sets of icosahedra may be arranged to cancel the remaining non-vanishing 
low-/ contributions. Our hope is that the energy bounds derived from this 
constrained problem provide a reasonable picture of the full unconstrained 
model. 

If there are icosahedra of fives, and icosahedra of sevens, the topo- 
logical constraint Eq. becomes 



- = 1 . (44) 
For a given configuration the energy is given by 

E""^ = Cni^R' + 12(2n+ - 1) E^r, (45) 

where the C coefficient is a function of 3(n^ + — 1) = 6n^ variables. Let 
us first choose a distinguished icosahedron with explicit coordinates 

{27rfc TT 2'7Tk 1 

(0, any), (7, -^)o<fc<4, (tt - 7, ^ + ^)o<fc<4, (tt, any) \ , (46) 

where 7 = cos~^(l/v^)- Each of the remaining icosahedra may then be 
parametrized by the set of three Euler angles necessary to bring them to the 
position described by Eq. 
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Figure 3: Particle configurations near a finite grain boundary. Circles repre- 
sent five-coordinated sites and squares represent seven-coordinated sites. 



The problem is to minimize over this set of Euler angles. We perform this 
minimization using a direction set algorithm [32|. From the results shown in 
Table |], it is clear that the energy coefficient C is reduced by the addition of 
defects. It is therefore favorable to form defects for a sufficiently small core 
energy. 



T 


Total 


CnIt 


a/R 


M 


1 


12 


0.60 


7 


12 


2 


36 


0.45 


0.09 


1791 


3 


60 


0.38 


0.06 


4031 


4 


84 


0.34 


0.03 


16124 


5 


108 


0.30 


0.02 


36279 


6 


132 


0.257 


0.02 


36279 



Table 2: Table of results for the minimum energy coefficient, as defined in 
Eq. (|i3|), obtained within the icosahedral approximation as a function of the 
number of icosahedral clusters of positive charge. The The penultimate 
column gives the average geodesic distance between neighboring charges. The 
last column gives the corresponding total number of particles, as estimated 
from Eq.dl). 



Another important issue is the precise arrangement of defects in the 
ground state. For n+ < 5 we find the remarkable appearance of finite grain 
boundaries - finite strings of interlaced fives and sevens, as depicted schemat- 
ically in Fig. ^. These grain boundaries are not always perfectly linear, al- 
though one does find alternating disclination chains clustered along geodesic 
line segments. Occasionally one finds dislocations, i.e. disclination pairs, 
displaced from this geodesic by a few lattice spacings. The ground state we 
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find, for tlie case n+ = 4, is depicted in Fig. ^ and nicely illustrates the 
above features. Note that although the local structure of these grain bound- 
aries mimics that expected for fiat space, the curvature of the sphere allows 
these linear structures to terminate, consistent with the discussion in the 
Introduction. 

For the case ?t,_|_ = 5 we still observe finite grain boundaries, but they 
show a tendency to spiral. This tendency becomes more pronounced for the 
case n+ = 6, where the finite strings evolve into more complicated structures. 

The picture emerging then within the icosahedral approximation is that 
small core energies favor a proliferation of defects. Below a critical number of 
defects of order 100, the ground state is well described by twelve finite grain 
boundaries, each one seeded by a defect in the original icosahedron. Above 
the critical number of defects the finite grain boundaries tend to branch and 
develop tentacles: the linear character of the pattern is lost. 

5.2 String dominated regime 

In this section we examine the relative orientation of the finite grain bound- 
aries discussed in the previous section. One ansatz is provided by a solution 
having the form depicted in Fig. ^ There is an axis of three-fold rotational 
(C3) symmetry at the center of the triangle formed by the geodesies connect- 
ing the three nearest-neighbor disclinations of the icosahedral array which 
forms the starting point of this variational ansatz. Finite grain boundaries 
are constructed by adding defects along the geodesic which joins the purely 
icosahedral sites with the center of Ca-symmetry. The midpoints of the grain 
boundaries form an icosahedron. The only free parameter in the model is 
the lattice spacing. This parameter may be fixed by minimizing the energy 
with respect to the lattice spacing, 

^ = 0. (47) 

da 

The interpretation of this extremal lattice spacing is discussed in Sec. ^ 

The results from this minimization are shown in Table ^ When the total 
number of defects is less than a critical value (approximately 110) this C3 
solution has energies slightly lower than those found within the icosahedral 
approximation. This is remarkable if one recalls that this Cs-symmetric 
solution is obtained by minimizing with respect to only one parameter, the 
lattice spacing. The results obtained from the icosahedral approximation 
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Figure 4: Six views of the ground state configuration for the icosahedral 
solution with seven sets of icosahedral defect clusters. The top figure in each 
column shows the north and south pole respectively. The subsequent views 
are obtained by successive rigid body 120° rotations of the entire sphere, 
using the right-hand rule, about an axis running from the north pole to the 
south pole. 
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Figure 5: Six views of the ground state configuration for the C3 solution with 
a large number of defects. The views are related as in Fig. ^. 
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Tntfll 

± KJ LfXl 


n 


n 
ti 


12 


60 

\J •\J\J 


'~V 

I 


36 


0.44 


0.121 


60 


0.37 


0.085 


84 


0.34 


0.062 


108 


0.32 


0.051 


132 


0.31 


0.042 


252 


0.28 


0.024 


492 


0.26 


0.012 


972 


0.255 


0.006 



Table 3: The minimum energy coefficient (see Eq. (0)) for the C3 solution, 
as a function of the total number of defects. The last column gives the lattice 
spacing a as determined from Eq. (^). 

itself are in rough agreement with this very simple ansatz, as apparent from 
Fig. p. Table |^ also makes clear that there is little gain in energy when the 
total number of defects exceeds the critical value, even in the limit of a very 
large number of defects. This is consistent with the picture that purely linear 
finite grain boundaries are replaced by more complicated structures when the 
number of defects is large. 

A more sophisticated treatment of strings, motivated by the discussion 
of dislocations in the Introduction, would build the grain boundaries from 
disclination dipoles with fixed size and then allow a variable spacing between 
these dislocations. We hope to pursue this approach in a future publication. 

5.3 Large number of defects 

We now have conclusive evidence that additional defects can lower the total 
energy of the system for small core energies. Defects will then proliferate and 
form highly complicated patterns. A detailed investigation of this regime is in 
progress, with complete results to be presented elsewhere. In this section we 
present one example of a branched structure that has lower energy than any 
of the linear structures considered so far. The structure we analyze consists 
of defects arranged in star patterns, or pentagonal buttons in the terminology 
of Toomre|T2| (see Fig. 0). To study these structures we construct rings of 
five-disclinations forming a pentagon with its center at the position of the 
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Figure 6: Six views of the ground state configuration for the superposition 
of the C3 solution and the icosahedral solution Fig. ^. The views are related 
as in Fig. ^. 
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icosahedron, as shown in Fig. |^. As free parameters we leave the radius of 
each ring, and the angle each ring forms with a given geodesic joining the 
center of a star to that of its neighbor. The topological constraint Eq. (p^ 
implies that there must be the same number of rings of fives as sevens. 



Figure 7: Example of a simple star defect with two rings, one of fives (circles) 
and one of sevens (squares). 

From Table ^ we see the energy for 132 defects is marginally lower than 
the corresponding value for the icosahedral approximation. For more defects 
the star clusters have significantly lower energy than the C3 solution. It is 
remarkable that all the disclinations in this ground state solution, other than 
the twelve seed disclinations, bind to form radial dislocations as illustrated in 
Fig. Furthermore, the relative orientation of the different rings conform 
to a rhombic tiling of the sphere consisting of 30 completely regular diamonds 
(the rhombic tricontahedron) , as shown in the bottom left picture of Fig. p. 
Note that by minimizing the defect elastic energy we obtain a dynamically 
generated particle spacing, for a fixed sphere radius, which optimizes the 
given structure. Further investigation of these pentagonal buttons, as well 
as other more involved structures, will be presented in the future. 



Total C a/R 



132 


0.255 


0.04 


252 


0.170 


0.025 



Table 4: Table of results for the minimum energy coefficient, as defined in 
Eq. (^) obtained within star defects as a function of the total number of 
defects. The last column gives the value of the particle spacing a. 

It is natural, at this point, to ponder the nature of the ground state 
in the limit of an infinite number of defects with vanishing core energy or, 
equivalently, in the limit R/a ^ 00. We have, in fact, already addressed this 
question in subsect. ^]^, where we proved that the only zero energy solution 
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Figure 8: Six views of the ground state configuration for four rings of pentag- 
onal buttons. The views are related as in Fig. |^. The bottom left view shows 
the associated rhombic tiling (the rhombic tricontahedron) of the sphere. 
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is an arrangement of defects {qi, {Oi, 4>i)}i=i ... satisfying Eq. (PT|). In this 
case the defect density would be fully rotational invariant and screen out the 
Gaussian curvature completely. Currently we find solutions that seem to be 
converging to this limiting case, but it is open as to how accurately one can 
achieve the desired limit C = 0. 

6 Instabilities of Icosahedral Lattices 

Our discussion so far has focused entirely on analyzing the distribution of 
topological defects on the sphere. We turn our attention now to the impli- 
cations for the underlying lattice structure. We thus take into account the 
regular six-fold coordinated nodes as well as the defects and examine the 
resultant lattices. 




Figure 9: The construction of a type (n,m) icosadeltahedral lattice. The 
filled circles indicate two nearest-neighbor five-fold disclinations. 

In the limit of large core energies our model predicts twelve disclinations 
forming an icosahedron. Lattices whose only defects are twelve positive discli- 
nations sitting at the vertices of an icosahedron may be constructed easily, 
since they are characterized by the path between two nearest-neighbor discli- 
nations. For a type [n, m) lattice this path consists of n straight steps from 
a given disclination, a 120° turn, and then m more straight steps to the 
nearest-neighbor disclination (see Fig. ^). The total number of particles M 
within this (n, m) (icosadeltahedral) lattice is 

M = 10(m^ + + mn) + 2 . (48) 

Within our model, the energy for these configurations has been computed in 
Eq. (PI). Since the core energy is sensitive to the short-distances properties 
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of the model, different icosadeltahedral lattices will have different energies, 
even for an arbitrarily large number of particles. 

Since most of the triangles in an icosadeltahedral lattice cannot be equi- 
lateral there is no uniquely defined lattice spacing. An average lattice spac- 
ing a can, however, be estimated. On the sphere, the distance between two 
nearest-neighbor disclinations is given by Rj (7 = cos~^(l/v^))- For the 
{n, 0)-case we have the relation Rj = n x a, and therefore 

a = ^ . (49) 
n 

Any other sensible way of estimating the lattice spacing, such as the size of a 
disclination dipole, should give a value of the same order. In the following we 
restrict ourselves to (n, 0) lattices for simplicity, but it is easy to generalize 
the formulas to arbitrary {n, m) icosadeltahedral lattices. 



0.3 ' — ■ — • — ■ — • — ■ — • — 

0.05 0.1 0.15 

a 



Figure 10: The C coefficient as a function of the lattice spacing for the C3- 
solution corresponding to a 3 length finite grain boundary. The straight line 
corresponds to the C-coefficient for a pure icosahedron. 

From previous sections we know that there will be a range of core energies 
for which the icosahedral lattices will be unstable to the formation of defects. 



To visualize this more clearly, let us take the Cs-solution of section q]2| for 
the case of finite grain boundaries having just three defects (a 5-7-5 con- 
figuration), and plot the C-coefiicient as a function of the lattice spacing. 
The result is shown in Fig. For sufficiently large lattice spacings the C 
coefficient exceeds that of a single icosahedral lattice. As the lattice spacing 
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is reduced, there is a critical particle spacing a*, such that the C coefficient 
of the C3 solution becomes smaller than that of the pure icosahedral lattice. 
Restricting ourselves to the (n, 0) icosadeltahedral configurations, we find 
from Fig. Eq. (H) and Eq. (|HD 

ajR ~ 0.2 ^ n(a,) = 5.5 ^ M(a,) = 305 , (50) 

where we exphcitly display the dependence on a^. This result implies that 
the pure (n, 0) icosadeltahedral lattice is unstable to the formation of defects 
for sufficiently small core energies and more than 316 particles. Alternatively, 
Eq. ( ^31) gives M = 363, consistent with the estimate above. Let us point out 
that the minimum of the energy occurs at lattice spacing = 0.121 (number 
of particles = 1256). For lattice spacings below Oc the disclinations 
will prefer to remain separated by a distance ac. This is accomplished by 
stringing six-fold coordinated particles between defects. Since lower energy 
configurations may be formed by allowing interpolating dislocations instead, 
we regard Mc as the maximum number of particles for which this particular 
structure is stable. The particle numbers quoted in the tables of Sec. |^ should 
be interpreted as the corresponding Mc. 

This shows that the minimum of the energy is attained for lattice spacings 
smaller than those necessary for stability. Finally, in the limit of vanishing 
lattice spacing, the Ca-solution becomes equivalent to a pure icosahedron, 
and both C-coefficients merge. 



7 The Thomson problem 

The Thomson problem can be stated as finding the ground state of an arbi- 
trary number of positive charges interacting through the usual 3d Coulomb 
potential, but with the further constraint that these charges must lie on a 
sphere. Since this problem falls within the universality class of our model it 
serves as a good testing ground. 

The Thomson problem has proven to be extremely difficult to analyze 
numerically, basically because of the large number of metastable states. Early 
analyses |TD[ showed that the ground state of the system for small numbers 



of charges was an icosadeltahedral lattice. Since a small number of charges 
corresponds to a large particle spacing this follows from our model as well. 
Some rules were also conjectured to decide on the true ground state when 
several icosadeltahedral lattices were possible for a given number of charges. 
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These rules could be examined within our model provided we are able to 
isolate the dependence of the core energy on the lattice type. 

Subsequent numerical work |Tl[] provided convincing evidence that the 
ground state, for a sufficiently large number of charges, does not have icosa- 
hedral symmetry. The critical number of charges for which additional defects 
arise seems to be around 400 |12], which is in agreement with our results. It 
is also found that these additional defects first arrange themselves into finite 
grain boundaries ||T^, as seen in our model. For more charges the ground 



state in the Thomson problem becomes very complex and the true ground 
state is not known. New configurations (one of the simplest being pentagonal 
buttons) appear to be energetically favorable in the early stages of this limit. 
This observation |]T2| is in agreement with our model as well. 

In the work of 0] (see also ||13|) it is also observed that the ground state 
energy for a large number of charges seems to converge to the energy that one 
would obtain in the unrealizable situation that all the charges are located on 
equilateral triangles. This limit corresponds to the defect density completely 
screening out the Gaussian curvature. The defect density therefore satisfies 
Eq. (|3TD , which we proved is the absolute minimum of our model in the limit 
of vanishing core energy. 

We think that the comparison of our model with the Thomson problem 
is very promising, but requires more detailed investigation to be addressed 
in the future. 



8 Discussion and Conclusions 

In this paper our first task was to propose and study an effective free en- 
ergy for disclination defects in particle arrays constrained to move on the 
surface of a two-dimensional sphere. The finite-temperature problem does 
not seem to be analytically solvable but we propose a discretized Laplacian 
Sine-Gordon model amenable to direct numerical methods. The structure 
of the ground state may, however, be studied analytically. This structure 
depends on the ratio of disclination core energies to the Young's modulus. 
On the sphere topology demands there be a total excess disclinicity charge 
of twelve. This excess charge can seed new ground state structures, com- 
pared to fiat space. For large core energies (or R/a < 36Ecore/ i'^^Koo^)) the 
disclinations arrange themselves to form an icosahedron. For intermediate 
core energies (i.e. R/a > 36 E core / i'^Koa'^)) grain boundaries develop which 
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terminate freely within the medium. The regime of still lower core energies, 
corresponding to R/a — oo, was found to be surprisingly complex - new 
defect arrangements make their appearance. 

Currently we are actively investigating the regime of small or vanish- 
ing defect core energy, including a detailed comparison of the predictions 
of our model with numerical results from the Thomson problem. A rigor- 
ous determination of the ground state for the Thomson problem is presently 
computationally prohibitive when the particle numbers exceed (9(500). Our 
methods enable us to reach particle numbers of 0(10, 000) or more with the 
same computational effort. 

Finally we believe that the rich symmetry structure underlying Eq. (|3lD 
may provide a direct analytic determination of the exact ground state in 
the limit of a large number of particles and further work in this direction is 
certainly warranted. 
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A The infinite radius limit 



In this Appendix we discuss in more detail the large R limit of the energy 
function Eq. (p9|) . It is readily seen from dimensional analysis and linearity 
that the x function of Eq. (^) scales like R^. 

Now consider a single isolated disclination qi located at point Pi on the 
sphere, together with a single dislocation, with Burgers vector b2, located at 
point P2- From Eq. ( ^9|) the total energy is 

i^ = ^(g?i?' + |b2|gi/(Pi,P2)i?+|b2riog(^)) ' 

where f{Pi,P2) is a function whose explicit form does not matter for the 
present analysis. The quadratic R dependence comes from the isolated discli- 
nation, the linear R dependence comes from the dislocation-disclination in- 
teraction and the logarithmic term comes from the dislocation energy. In the 
infinite-radius limit of the sphere we see, therefore, that the various defect 
energies scale identically to those in a flat space system of size R [^, |2B| . 
The nature of the ground state, however, is dramatically changed. 



B The biharmonic operator on the sphere 

The evaluation of the inverse biharmonic operator on the sphere is rather 
tedious. We outline the steps for a sphere of unit radius. The simplest 
approach is to first compute the inverse harmonic operator. The sum over m 
is performed as in Eq. (|26|). The result involves only Legendre polynomials 
and is 

The sums over / may be performed using the identities 

(1 - 2ua; + 



1=0 

and 



1 1 
1 Piix) = / du- KmK (53) 
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The resultant integrals are readily evaluated, yielding 

r(x) = -log(i^)-l. (55) 
The inverse biharmonic operator now follows from the result Ax{x) = r{x). 



C The evaluation of % 

Performing a trivial integration by parts in Eq. (|28|) , we get 

xiP) = 1 + t(log(t) - 1) + r dz^^ , (56) 

Jo I — z 

with t = l^l^a^. The last integral may be expressed via a change of variables 

as 

z logz xe~ 



X 



dz- = - / dx— — - , (57) 

I — z Ju — 1 

with u = ln{l/t). Finally one can expand for small u and large u as 



n+1' n+l 



where Bk are the Bernoulli numbers. These expansions are very useful as 
they allow a numerical evaluation of x with arbitrary precision with negligible 
computational time. In fact we save, on average, a factor of 2000 in time 
compared to a direct evaluation of the integral using Romberg integration. 
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